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Abstract —Initial ranging constitutes a part of the synchroniza¬ 
tion procedure employed hy the wireless communication stan¬ 
dards. This allows the base station (BS) to detect the suhscriher 
stations (SS) that are willing to commence communication. In 
addition, the ranging process allows the BS to estimate the 
uplink channel parameters of these SSs. Accurate estimation 
of these parameters are crucial as they ensure that the uplink 
signals from all the SSs arrive at the BS synchronously and 
approximately at the same power level. However, this detection 
and estimation problem turns out to be very challenging when 
multiple users initiate the ranging procedure at the same time. 
We address this issue by exploiting the underlying sparsity 
of the estimation problem. We propose a fast sparse signal 
recovery approach to improve the ranging performance in multi¬ 
user environment. Compared to the standard correlation based 
techniques, our method shows a clear Improvement in ranging 
code detection, timing offset and channel power estimation. 
Although this method has been developed around the WiMAX 
standard, the underlying principles apply to other OFDM based 
standards as well. 

Index Terms —Initial ranging, code detection, sparse represen¬ 
tation, OFDMA. 


I. Introduction 

The orthogonal frequency-division multiple access 
(OFDMA) scheme has been adopted by the IEEE WiMAX 
standard H). To maintain orthogonality among the subcarriers 
of different users, and to avoid the occurrence of multiple- 
access interference (MAI), the uplink signals arriving at 
the BS must be aligned to the local time and frequency 
references. For this purpose, the IEEE WiMAX standard 
enforces a network entry procedure called initial ranging 
(IR). The ranging process starts with the allocation of a 
pre-defined set of subcarriers by the BS in specific time slots, 
which is known as a ranging opportunity. The SSs which 
wish to commence communication with the BS, referred to 
as the ranging terminals (RTs), can take this opportunity 
by modulating a randomly selected ranging code onto the 
allocated subcarriers. Due to different position of RTs within 
the radio coverage area, ranging signals transmitted by 
different RTs arrive at the BS with their specific transmission 
time delay. At the receiving side, the BS is required to extract 
timing and power information for each detected code and 
inform the RTs about the extracted information. 

Multiuser code separation as well as their timing and power 
estimation are the main tasks of the IR process ||2|. The corre¬ 
lation based approach proposed in lO, is based on the principle 
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that a time delay can be represented by a phase shift in the 
frequency domain. Lee replaces the WiMAX ranging codes 
by a set of generalized chirp-like polyphase sequences to get a 
more accurate timing estimate. The work in fS) demonstrates 
that the frequency-domain correlation approach outperforms 
its time domain counterpart. However, the methods a-a 
simply treat the MAI as a noise, which results in performance 
degradation in multiuser environment. An approach different 
from the IEEE 802.16 standards has been proposed in Q 
to overcome the MAI problem in code detection. The idea 
is to allocate a small number of subcarriers to each ranging 
opportunity so that most of the RTs are expected to transmit on 
disjoint sets of subcarriers resulting a minimum level of MAI. 
However, the reduction of the number of effective subcarrier 
for each user results in the degradation of timing estimation 
performance 121 . A similar approach has been proposed in ||7l 
for channel synchronization. This method assumes that the 
uplink signals are transmitted over disjoint subcarriers, and 
the receivers use filter banks to separate multiuser codes. The 
work in fSl improves the ranging performance by dividing 
the ranging signals into several groups with each group 
being transmitted over exclusively assigned subcarriers. The 
concept of successive interference cancellation (SIC) has been 
employed in lH, 13. The method proposed in 121 works in an 
iterative fashion where the strongest path of each active RT 
is detected and is removed from the received signal and the 
resulting signal is used in succeeding iterations. The work l9l 
differs from m in the aspect is that in l9l one additional 
ranging signal is detected at each iteration instead of a single 
multipath component. 

Sparse signal representation Qol, O has many potential 
applications, including spectral analysis ina, channel estima¬ 
tion 03 etc. In this paper, we pose the problem of user 
detection and channel estimation in ranging process as a 
sparse recovery problem. This approach is founded on the two 
observations: 

■ The number of ranging subscriber stations is much 
smaller than the number of ranging codes. 

■ Only a few channel taps are of some noticeable magni¬ 
tude, and the remaining vast majority of channel taps are 
of negligible magnitude. 

We first develop a signal model that allows us to exploit 
the above facts and pose the ranging problem as a sparse 
recovery problem. By applying the standard sparse recover 
methods to solve the sparse recovery problem arising in 
ranging problem, we found that associated computation time is 
too high compared to what is typically allowed by the WiMAX 
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Standard. For this reason, we propose a new algorithm to 
solve the sparse-ranging problem. The proposed algorithm 
combines two different types of sparse recovery algorithms to 
provide a time efficient solution. In particular, we apply a non- 
convex sparse recovery algorithm which has fast convergence 
property if it initialized sufficiently close to the final solution. 
To obtain a good initial input, we consider another convex 
optimization algorithm that can provide a good estimate of 
the final solution after a few iterations. The proposed handover 
algorithm computes most of its multiplications by using Fast 
Fourier Transform (FFT). Due to noise contribution in the 
ranging signal, the recovered signal from the algorithm may 
not be truly sparse. For this reason, we analyse the statistical 
property of the recovery error. This analysis is used to formu¬ 
late a systematic hypothesis test to detect the codes and the 
associated timing offset. 

Notation: Lowercase boldface letters denote vectors. The 
A:-th component of a vector x is denoted by [a;]^. Uppercase 
boldface letters are used to represent matrices. For a matrix 
A, we use [A]y to denote the element of A at its i-th row 
and j-th column. The Z denotes the set of all integers. We use 
(A)*, and (A)''' to denote complex conjugate transpose and 
transpose of the matrix respectively. and !„ are column 
vectors of size n with all O’s and I’s respectively. 

II. Signal Model 
A. Single ranging terminal 

Consider an uplink OFDMA system with N subcanders. 
Thus, each OFDM symbol contains N data symbols. In 
addition, an OFDM symbol must contain a cyclic prehx. Let 
the length of the cyclic prehx be Ng. Therefore, the length of 
a OFDM symbol is 

N = N + Ng. 

Now consider a particular transmitter T in this system. 
Let u{k), /c S Z be the sequence of the channel symbols 
transmitted by T. These channel symbols are grouped in 
OFDM symbols. Without any loss of generality we assume 
that the nth OFDM symbol consists of u{k), k = 
nN ,nN +1,... ,nN + N —l.The n-th OFDM symbol vector 
is written compactly as 

Un = [u{nN) u{nN + l) ■■■ u{nN + N — 1) Y. (1) 

Let us denote the impulse response coefficients between the 
transmitter T and the base station as h{p),p = 0,1,..., P — 1. 
Consequently, the contribution v{k), fc S Z of transmitter T in 
the signal received by the base station is given by 

AT-l 

v{k) = hpu{k — p — d). (2) 

p—0 

In practice, the value of N is more than 200, and typically 
h{p) = 0 for all p > 50. The delay d depends on the distance 
between T and BS. One purpose of the initial ranging process 
is to estimate d so that the transmitter can align its transmission 
with the frame boundaries of the base station. In a cellular 
communication architecture d cannot be arbitrarily large. It 


must be bounded. In IEEE802.16 the cell radius is chosen 
such that d < N. 

IEEE802.16 identifies some specific uplink subcarriers as 
‘ranging subchannels’. In the sequel we assume that there are 
M subcarriers in the group of ranging subchannels, and denote 
their indices by {jm ■ rn — 1,2,...,M}. The value of M 
in IEEE802.16 is 144. In IEEE802.16 a ‘ranging opportunity’ 
spans over two consecutive OEDM symbols when the ‘ranging 
terminals’, who wish start communicating via the base station, 
can send their ranging codes. The ranging codes must be sent 
via the ‘ranging subchannels’. 

Consider a ranging opportunity consisting of the OEDM 
symbols n — 1 and n, respectively. Suppose T is a ranging 
terminal who wants to use this ranging opportunity. According 
to IEEE802.16, T must construct the OEDM symbols Un-i 
and Un as follows. Eirst it chooses a column q of a pre- 
specihed M x G code matrix 

C = [ Cl C2 • • • Cc] 

uniformly at random. Thus the probability that f is a given 
integer in the set {1, 2,..., G} is 1/G. IEEE802.16 specifica¬ 
tion defines C preciely. Hence it is known to the BS and all 
the transmitters including T. After choosing i the terminal T 
calculates the numbers 

1 ^ 

Sg = -j= ^ [C]m.i exp{i2Trjmq/N}, g = 0,1,..., A - 1, 

^ m—1 

(3) 

where [C\m^i denotes the m th component of ci, which is 
also the element at the m th row and £ th column of C. The 
operation (O can be seen as the process of modulating the 
jm th subcaiTier by Cm,i 7 and modulating the non-ranging 
subcarriers with 0. In practice this computation is carried out 
using the lEFT algorithm, and is compactly given as 

s := [ So Si ••• SAT-i ]t = F*©tcf, (4) 

where F is the N x N FFT matrix such that 

[F]fem = exp{-27ri(fc - l)(m - 1)/A}/v/A, i = \/^. 

The matrix © is an M x A row selector matrix such that the 
m th row of © is the jm th row of the N x N identity matrix. 
Recall that jm,rn = 1,2, ...,M are the indices associated 
with the ranging subchannels. The process of modulating the 
subcarriers via IFFT is equivalent of pre-multiplication by 
F* in (IHi. The premultipaction of q by ©t implies that 
only the ranging subcanders are modulated by the appropriate 
components of the ranging code cg. 

Using s the transmitter constructs the OFDM symbols m„_i 
and Un as 

Un-l = [ SN-Ng ■■■ SN-2 SN-1 So Si ••• SjV-l]"'’, 

Mra = [ So Si ••• SN-1 So Si ••• SAr^-i]'''. (5) 

To detect the ranging codes the base station works with the 
first A samples of the nth received OFDM symbol. According 
to the standard practice in OFDMA, the BS computes an FFT 
of these samples and then examines the data received in the 
ranging subchannels. 
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We first find the contribution 


Vn = [v{nN) v{nN + l) ■■■ r!(n-/V + — 1) ]''' 


of T in the first N samples of the nth OFDM symbol received 
by the base station in terms of s and 


h = [ ho hi 


hN-i ]^- 


( 6 ) 


Using O and 0 note that both m„_i and are linear 
functions of s. In particular, 


u{nN — k) 


s{k-N) -N <k<-N, 
s{k) —N < k <0, 

s{N-k), 0<k<N, 

s(27V-/c), N<k<N. 


In the following we denote the circular shift operator by (•). 
For instance, the circularly shifted version of s by /c places is 
given as 


^iik)'=[sN-k ■■■ SN-1 So Si ■■■ SAr_fe_i (8) 

A few steps of algebra using (|2]i and 0, and using the fact 
that d < N, we get 


= (9) 

where H is an iV x cyclic matrix 

H = [ ^4,(0) ^4,(1) ^4,(iv-i) ]• (10) 

Expressions (|9]l and ( fTOl l involving circularly shifted versions 
of s and h are typical of OFDM, and are resulted from the 
way the OFDM symbols and m„_i are constructed in 0. 
This special construction allows us to exploit the identity na 

Fa;|(fc) = diag(F(:, k + l))Fa; 

= diag(Ftc)F(:,fc +1), fc = 0,1,..., A^ - 1, (11) 

satisfied by the FFT matrix F, and any vector x. Note that 
we use the standard MATLAB notation F(:,fc) to denote the 
fcth column of F, and diag(a;) to denote the diagonal matrix 
such that [diag(a;)]fc fe = Recall that for detecting the 
ranging codes the BS must compute an FFT of t;„, and then 
extract the data received in the ranging subchannels. The 
FFT of Vn is Fit„. Hence, the data received from T in the 
ranging subchannels is given by premultiplying Frt„ by the 
row selector matrix 0. Using (fTOl i and (fTTT i it follows that 

FH = diag(F/r)F. 

Note that h := Fh is the FFT of h. Hence by 0 and (fTTT i if 
follows that 

Ft;„ = FHs|(d) = diag(/i) Fs^(^d) 

= diag(/r) diag(F(:, d + 1)) Fs 
= diag(Fa) diag(F(:, d + 1)) h 
= diag(Fs)Fh.4,(d). (12) 

It is well-known that F is an unitary matrix, i.e., F*F = 
FF* = I. Hence 0 implies Fs = Hence the 

data received by the BS at the ranging subchannels due to 
transmission of T is given by 

0Ft;„ = E|h,4,(d), Ef = 0 diag(0TQ) F. (13) 


Each of the matrices E^,£ = 1, 2,..., G is of size M x N, 
and is known because Ci is known. 

Typically, we know a number P such that \hk\ = 0 
for k > P. At this point we emphasize that OEDM can 
effectively equalize the inter-symbol interference effects only 
when P < Ng. Thus the existence of the upper bound P is a 
key assumption in OEDM. In addition, the cell radius gives an 
upper bound D on d. Hence d+P < D + Ng. By construction 
of ft 4 .(d), we know only first 79 -f P of its rows are non-zero. 
Hence it is fine to truncate to a D + Ng dimensional 

vector, and thus it is enough to work with only first D + Ng 
columns of E^. 


B. Multiple ranging terminals 

So far we have considered a single ranging terminal, and 
in (foi l we have quantified its contribution to the data vector 
received by the BS in the ranging subchannels. In this section, 
we generalize the analysis for multiple ranging terminals, and 
account of the receiver noise at the base station. 

Suppose that the code q is chosen and transmitted by iV^ 
number of ranging terminals. This means that the total number 
ranging terminals in the system is 7T = A^i -I- A ^2 + • ■ ■ + Nq. 
We emphasize that Ng is a random quantity for a given £. 
Typically 7T < 10, G = 256, and the probability that Ng > 1 
is 

Prob{iV^ > 1} = 1 - (1 - 1/G)^ - K{1 - l/G)^-i(l/G), 

which is a very small number. Hence the probability that two 
or more ranging terminals will collide by selecting the same 
code is very small. However for the sake of generality we do 
not exclude that possibility. 

Let h\ ' and dg^k-, k — 1,2,Ng denote the channel 
impulse response vector and the delay of the fcth ranging 
terminal transmitting the code cg. Then by the principle of 
superposition the data vector y received by the BS at the 
ranging subchannels is given by, see (fOl i 

G 

y ^'^Fghg + e, ( 14 ) 

i=i 

where e is the additive receiver noise, and 

hg = l '^kLlVA (15) 

\ 0, Ng = 0. 

is the combined channel vector for all the ranging terminals 
transmitting the code cg. Note that the power received by the 
BS corresponding to the ranging code cg is given by El: 

Ff = h*ghg. 


HI. Estimation of ranging information 
A. Formal problem Statement 

Given y the signal model in (fT4l i the BS needs to 

1) Find the set C = {i : Tg ^ 0}; 

2) For every £ G C find F^ and dg^i assuming Ng = 1, \/£ G 
C. 
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In some rare cases > 1 for some i. In this case the 
ranging requests of the users who chose ci simultaneously 
would collide. Nevertheless, as we see later, the BS would 
detect that cg was transmitted among others. A IEEE802.16 
base station allocates some small bandwidth corresponding 
to every detected ranging code. The ranging terminals use 
this ‘grant’ to transmit their buffer status information, and 
expect to obtain some adequate amount of bandwidth from the 
BS to commence data transmission. However, if two ranging 
terminals, say T1 and T2, chose the same code cg for ranging 
request, then the estimate of E^ and dg^i obtained by the BS 
during the ranging process would have no physical meaning. In 
addition, the bandwidth request from T1 and T2 will collide 
again. Consequently, the BS will not be able to decode the 
bandwidth request data from T1 and T2. In such a scenario 
a IEEE802.16 BS does not allocate any further bandwidth 
corresponding to code cg, and after a timeout period T1 and 
T2 commence the ranging process again ||2l, ||3l. 

Recall that the first d components of are zero, see 

dSll. Hence by construction of hg in (flSl l. the index of the first 
nonzero component of is 1 + dg, where 

dg = min dg k, i € C. 
k(g{l,2,...,Ne} 

Clearly, if Ng = 1, then dg = dg^g. Eor this reason we 
propose to estimate dg as the timing offset corresponding to 
an f S £. When Ng = 1 this estimate is consistent with our 
requirements. On the other hand, if Ng > 1, this estimate will 
have no practical relevance, for, as discussed above, BS will 
reject cg in the bandwidth request stage. 

B. Sparse recovery framework 

Recall that for any f only first 

Ni-.= D^ Ng (16) 

components of hg are non-zero. Hence 

^ghg = Y.g{:,Ng) hg{l:Ng). 

Note that we use Matlab notation E^(:, 1 ; Nf) to denote 
the submatrix of Eg formed by taking its hrst columns. 
Similarly, hg{l : Ng) denotes the vector formed by taking the 
first Ng components of hg. Then we can write (fl4l l as 

y = Ax -f e, (17) 

where 

X[ h-[{l : Ng) hl{l:Ng) ••• ] t , 

A=[Eg{:,l.Ng) E2i:,l-.Ng) ■■■ Eg(:, 1 : A^i) ]. 

Note that by definition of Eg in (fOT l. A is a known matrix. 
On the other hand x and e are unknowns. Typically, the total 
number of ranging terminals K = Yl't=g implying 

Ng = 0 (and therefore hg = 0) for a vast majority of the values 
^ S {1, 2,..., G}. This makes x very sparse. This observation 
motivates a sparse recovery framework for solving the ranging 
problem. 

We propose to estimate a sparse vector x that is consistent 
with ([I7]i. There are many reliable algorithms for solving 


such sparse estimation problems Eoi, im. Denote the sparse 
estimate by x. Then the BS can extract the required ranging 
information as follows. Partition x into G number of sub¬ 
vectors: 

x = [hl hi ■■■ hoY, 

where each hg is of length Ng. Then we declare £ G C only 
if 11/if 11 f 0 and the index of the first nonzero component of 
hg leads to an estimate of dg. 

C. Background on sparse recovery methods 

If e = 0, then the ideal way to reconstruct a sparse x from 
y requires solving 

a;* = argmin ||i;||o subject to y = Av, (18) 

V 

where ||it||o, which denotes the fp norm of a vector v, is 
simply the number of non-zero components in v. Thus the idea 
is to find x^, with the smallest number of non-zero components 
satisfying y = Aa:*. The unique representation theorem ifTSl 
ensures that under mild technical conditions there is a unique 
a;* with ||x*||o < M/2 satisfying y = Aa:*. 

However, (fTSl l is combinatorial in nature M- The most 
popular alternative approach for relaxing (fTSl l is called Basis 
Pursuit (BP) iflTl . lIT^ . where the fp norm in (fT¥l l is replaced 
by .fi norm: 

x^, = argmin ||t;||i subject to y = Av. (19) 

V 

Here 

GNi 

ll^lli := IMfel- 

k^l 

BP can be posed as a linear program El over second order 
cones, and can be solved in polynomial time. In addition, it 
has been shown in El that BP recovers the sparsest solution 
to y = Ax with a very high probability. 

The above simple ideas can be adapted quite well even when 
e f 0 ifTOl . However, the existing algorithms for solving ET i 
are unable to converge to a satisfactory solution within the 
time-frame available to solve the ranging problem in practice. 
In the sequel we propose a new approach to overcome this 
hurdle. 

Our approach blends the nice properties of the £q and £g 
methods. The so called £q approximation methods 
are known to converge very fast if initialized sufficiently close 
to the final solution. But these methods being non-convex, may 
often get trapped in some local optimal point when initialized 
far away from the final solution. The £g methods being convex, 
does not have the local optima problem, but typically take a 
large number of iterations for convergence. Therefore, we aim 
to to start with an £g approach and then handover to an £q 
approach when the solution is ‘sufficiently close’. In particular, 
we introduce a new £g norm minimization method that can 
obtain a rough estimate of x in only a few iterations, and then 
handover to a reliable fo-^pproximation algorithm GOl, ED. 
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D. l\ optimization algorithm 

In this paper we propose to solve a special dual of ( fT9l l. 
This dual formulation relies on the theory of minimum norm 
problems in Banach spaces. Given element u G we 

define the infinite norm as 

||m||oo = max IMfcl- 

feG{l,2,...,GAfi} 

Let us define the bilinear from (•, •) as 

{v,u) := —{u*v + v*u). 

Note that this bilinear form maps x onto R. By 

Holder’s inequality it follows that 

< ||r;||i||M||oo, (20) 

provided that both ||rt||i and ||m||oo exist. In addition, when 
we have an equality in (l20l l. then we say v and u are aligned. 
The condition for alignment can be verified to be as follows 

ED: 

Proposition 1: Let 

K = { k : |[M]fe| = ||m||oo}- 

Then {v,u) = ||r;||i||tt||oo only if 

• = 0, for all fc ^ /C; 

• For all fc G /C it holds that = pfeConj([tt]fc)/|[M]fc| 
for some non-negative number p,k G M. 

We are now ready to state the key result allowing us to 
formulate a convenient dual of ( fTOl l. 

Theorem 1: Let us define the sets 

V={t;: At; = y}, U = {g : ||A*g|U < 1}. 

Then 

min||t’||i = max-(y*g-t-g*y). (21) 

z 

In addition, let tc* be the solution to the optimization problem 
in the left hand side of (ISlT i. and let g^, be the solution to the 
optimization problem in the right hand side of (ISTT i. Then 

(a;„A*g,) = ||a;,||i ||A*gJU. (22) 

Proof: See m 

Note that g is of significantly smaller size compared to x. 
Computationally it is a lot more economical to solve the 
dual in the right hand side of (ISTT i. and apply the alignment 
condition (l22l l to recover x* from g*. We write the dual 
problem as 

= argmax + y*g) (23) 

g 2 

subject to g*aia*g <1, i = 1,... GNi (24) 

where the i th column of A is denoted by ai. We wish to 
solve it via a primal-dual algorithm ll23l . Therefore we write 
the Lagrangian associated with (l2^ - (l24ll : 

z) = \{9*y + y*g - g*Adiag( 2 ;)A*g -f Vz} (25) 


where 1 is a GN^ dimensional vector of all ones, and 
z/2 is the real-valued vector of Largange multipliers. Each 
component of 2 ; must be non-negative, and we denote it by 
2 ; > 0. Now it is a standard result in the theory of least squares 
that 

argmaxL(g,z) = [A diag( 2 )A*]"^y, (26) 

maxL(g, z) = ^{Vz + y*[A diag( 2 ;)A*]"^y}. (27) 

Therefore, we can obtain 2 ;* by solving the Lagrangian dual 
of (|23ll-(|2l; 

2 :, = argmax l;{Vz + y*[A diag( 2 :)A*]“^y} 
z 2 

subject to 2 ; > 0. (28) 

In reality, a primal-dual algorithm would solve (l2^ - (l24l) and 
its Lagrangian dual (l28l l together by finding the solution to the 
Karush-Kuhn-Tucker (KKT) conditions 

y = A diag( 2 ;*) A*g^, (29) 

W^>0, i = f...,GNi, (30) 

(l-g:a,<gj >0, i = l,...,GAi, (31) 

[z,]^ il-g:a,a*gP = 0, t = l,...,GN,. (32) 

Equation ( |29] ) follows from (|26]). The inequalities (O and (1311) 
must hold because the constraints in (l24l i and (l28T l must hold. 
Equation (1321) is the complementary slackness condition which 
says that for any i either (1 — g%aia*gf} = 0 or [zf = 0. 
These relations can be used to calculate £c» from 2 :* and g„. 

Proposition 2: The optimal solution x, of (fTOl l is given in 
terms of z* and g, as 

X* = diag(z,) A*g,,. (33) 

Proof: Note that by setting x* as in (l3^ we do satisfy 
y = Ax*. It remains to verify the alignment condition (1221) . 

Now, it must hold that ||A*g*||oo = 1- This is because 
if ||A*g*||oo < 1, then we could always multiply g* by a 
suitable real valued scalar k > 1 such that ||A*(Kg^)||oo = 1, 
and 

aly + y*9* < {^gSy + y*{ngj, 

leading to a contradiction. 

Let us define 

IC = {i: gla,a*g^ = 1}. 

Note that K. is nonempty since ||A*g*||oo = 1. However, the 
complementary slackness condition (l32l l implies that \zf\i = 0 
for all i K.. Hence 

i£K, i^K, 

The last equality follows because by definition of 1C we have 
\a*g\ = 1 for all i G 1C. Then 

g*Ax* = ^g>ia-g[z*]i = = ||tc*||i 
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Hence we can verify that the alignment condition 


(a;,, 


+ a;*A*gJ = 1 x ||a;*||i 
||A*gJ|oo||a;*||i 


holds, and thereby the proof is complete. 


We wish to find a numerical method to solve the KKT equa¬ 
tions (I29b- (l32l) . which are nonlinear simultaneous equations in 
2 * and To solve the KKT equations using the primal-dual 
method, one can relax the complementary slackness condition 
in (l32l) to 


N*]* (1 - i = l,...,GNi. (35) 

where g > 0. The value of p decreases as we progress through 
the iterations of the primal-dual algorithm. The standard way 
to handle the modified KKT equations daiii-iEB and (El is 
to use the Newton’s approach. A derivation of the primal-dual 
algorithm can be found in ll2l . and in our case it reduces to 
the form summarized in Table-U where we define the function 


Mg) = g*a^a*g - 1, for i = l,2,---GNi (36) 

with f{g) = [fi{g) / 2 (g) • • • /GWi(g)]'^ and vectors g and b 
such that 


q = A*g 

= for * = l,...,GAi. 

Mg) 

The Jacobean matrix of /(g) turns out to be 


(37) 


J= [aiMg aaaag---aGWiaGATiS]^ = diag(g)*A*. 


(38) 


The algorithm in Table-H] terminates when a rough estimation 
of X has been obtained. We use the following procedure at 
every iteration to check whether the algorithm has yield a 
sparse enough estimate of x. Note that x is sparse and hence 
most of its energy should be concentrated on a few number 
of its components. In fact, it has been shown in ifTSl that any 
sparse recovery algorithm can perform well when the energy 
of X concentrates within M/2 number of its components. At 
Step-5 in Table-U we compute x. Let x be the thresholded 
vector constructed from x by retaining M/2 most significant 
components of x while setting others to zero. Now compute 



The value of k can be used as an indicator to measure the 
proximity of x to the actual x, while k ^ 1 indicates that 
X may be very close to x. However, the algorithm in TableU 
targets to produce only a rough estimate of x and hence we do 
not wait until k = 1. In particular, the termination value of k 
trades-off computational complexity with estimation accuracy. 
A larger terminating k brings x closer to x for increased 
number of iterations. The value of k is being an interesting 
design parameter. We have observed that the primal-dual 
algorithm can achieve a k > 0.6 in only 4 — 5 iterations. 
The reason behind taking small number of iterations may be 
due to the formulation of the primal direction search step i.e.. 


TABLE I 

Primal-dual algorithm 


Initialization 

1. Set g = 0, z = 0. and parameters /i, 0 < o < 1. 

repeat 

2. Compute primal-dual search directions: 

Ag = {A diag(z)A* — J*SJ)~^(—y -I- /rj*6) 

Az = —(z + gb + SJAg) 
where, S = diag(diag(/(g))-lz). 

3. Find 0 < s < 1 such that: 

3^- fiig + S'^ 9 ) < (3,z(*) + sAz{i) > 0; Vi. 

3b. The norm of residuals has decreased sufficiently: 
||T^(g -I- sAg,z + sAz )||2 < (1 - os)||v(g, z)|| 2 . 

4. Set g = g sAg, z = z sAz. 

5. Compute x = diag( 2 ;)A* 5 . 

6. Set fi = afi. 

until (A rough estimation of x has not been obtained) 


Ag (see Step 2, Table-U). The popular non-convex iterative re¬ 
weighted least square (IRLS) based algorithms like FOCUSS 
na, ISLO lIMI use similar formulations to update the estimate 
of actual sparse signal (x) in every iterate. It is well known that 
if IRLS algorithms can avoid local minima then they provide 
close estimate of actual signal in a small number of iterations. 


E. Smoothed fg minimization 

We use the rough estimate of x obtained by the primal- 
dual method to initialize the improved smoothed £q (ISLO) 
algorithm ll20ll . The ISLO algorithm is described below. Dehne 
the Gaussian functions, 

f^{a)=e~M. (39) 

Then it is readily verified ll20l that, as cr —)► 0, the function 

GNi 

Pcr{x) = fcr{[x]t) (40) 

t=l 

behaves like GNi — ||a;||o, motivating the following approxi¬ 
mate reformulation of (fTsT i: 

x^{a) := argmin — subject to g = Av, (41) 

V 

while taking cr —)■ 0. Like ||u||o, the function F^{v) has 
many local minima for a small a. Hence, one solves (EB 
for a large tr initially, and successively decrease cr using a 
small factor and solve (l4Tli repetitively. Finally solves (l4Tli for 
cr = (To, where ao is a small positive number. The work in 
II 20 II proposed a systematic way to choose erg. In presence of 
noise, i.e. when e 0 in (fTTl) . the following optimization for 
ISLO has been considered in mi: 

a:* (cr) = arg min L^{v), (42) 

V 

La{v).= -F„{v)+ ^\\y - Av\\l, 

where A > 0 depends on noise level. A Gauss-Newton type 
convex-concave procedure is used in ifT^ to minimize for 
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TABLE II 
ISLO Algorithm 


Input: (Tat 

Initialization 

1. Set fT = (Tst, A € [1,100] and £ [0,1), 

i = 0,0-0 6 [0.1,10-^]. 

repeat 

2. Set /3 = 1. 

3. while La{l3({x^^^) + (1 — /3)a:(*)} > La-(x^'^'>) 

13 = 7/3. 
end 

4. £c(*+i) = + (1 - Set ^ ^ + 1. 

5. If 1 12 < rju then a = pa. 

while a > aQ. 


a fixed a. A detailed description of convergence properties 
of dnii and (l42]) can be found in ifTSll . II 20 II . In particular, the 
following Lemma gives a direction to minimize for a fixed 
a. 

Lemma 1: ifT^ Define the mapping C, : 
such that 


C(*) =A[fL<,(*)/a2 + AA*A] ^ A*y, 
where Wa{x) is a diagonal matrix: 

/<t([*]i) ■■■ 0 


Wct{x) = 


f<T{[x]GNi) 


(43) 


(44) 


Then x*(ct) = C{a:*(tT)}. In addition, for any u there exits a 
real-valued scalar ^ > 0 such that 


L^{k({u) + {1 - ^)u} < L„{u). (45) 

Lemma [T] reveals the fact that x^a) = ^od 

motivates a fixed-point iteration approach to find x^a) by 
solving the equation x = C{®}- Furthermore, the Lemma 
provides a direction such that L^{x) is decreasing along 

The ISLO algorithm is given in Table-HIl Here denotes 
the value of x updated at the i th iteration. The procedure of 
choosing x^^'^ and Ust will be described in the next section. 
The value of A in (l42ll controls the distance between y and 
Aa;*(cr). A small value of A allows \\y — Aa;*(cr)||| to be 
larger. Note that according to (fTTl i. the term \\y — Aa:||| is 
equal to the power of measurement noise. Hence, we should 
choose a small value for A when noise variance is larger. A 
procedure for choosing the value of A is described in ll2^ . 
The popular choice of A = 0.1|| A*y||oo. For minimizing L^, 
we use C{x) — x as the descent direction. The 7 is a standard 
backtracking line search parameter ll23l . The inner-iteration 
for minimizing for a given a terminates when — 

< i?cr, (see Step 5). We then update a = pa. The 
work in ll20l describes a procedure for choosing the values 
of T] and p. In particular, it has been shown that the ISLO 
remains insensitive to the value of the parameters if we choose 
T] G [0.1,0.7] and p G [0.2, 0.9]. In this work, we set 77 = 0.5 
and p = 0.3. The stopping criterion of ISLO is based on a 


small value of a denoted by ctq which depends on the noise 
level. A wide range of numerical simulations in noisy cases 
(5 dB to 20 dB SNR) suggest that ao = 0.001 is good choice. 


E The Elandover algorithm 

The proposed handover algorithm starts with the £ 1 - 
optimization in Table-|I| However, we allow the algorithm to 
provide only a rough estimation of x. Once a rough estimate 
X is obtained it is used as initial x^^^ for ISLO (see Table- 
Hill. The value of a^t for ISLO can be found in the following 
way. Assume that x is the minimizer of L^-^^ (v) in (l42li . Then 
according to Lemmall] and (|4^ we have 


X = C{&} 

= A {W„{x)lalt + AA*A] K*y. (46) 


We need to solve (EUl for a St- However, the equality in 
(l46l l may not hold in practice. Then the value of a^^ can be 
approximated by: 


= argmin 


- AA*(y - Ax) 


(47) 


The optimization problem is non-convex, but one dimensional. 
An interior trust region algorithm ll25l has been applied to 
estimate ast from dlTll where the initial value of a is set to 
maxj(|[x]j|) (see ifT^ for justification). 

Suppose the final output obtained from ISLO is x. Due to 
noise contribution in (ITtI) . the estimate x may not exact copy 
of X, and hence x will have spurious peaks. As a result, x 
needs thresholding to perform the code detection. We shall 
develop a thresholding procedure in the Section jlll-HI Let 
X be the thresholded vector constructed from x. The BS 
can detect the active ranging codes and corresponding timing 
offsets from x by using the procedure described in Section- 
IIII-BI To obtain the estimate of channel impulse response 
(CIP), let us partition x = [h[ h\ ■ ■ ■ Iiq j^. The estimate 
of CIP corresponding to the £-th active code is hg. 


G. Computational Complexity Analysis 

In this section, we shall analyse computational complexity 
of the ISLO algorithm. A similar procedure can be followed to 
analyze the complexity of the £i algorithm in TableU As can 
be seen in LemmaU] the major fraction of the computation for 
ISLO is involved in computing C(x) in ( l4Tt . However using 
the matrix inversion lemma in (l43T l one can verify that 

C(i) = W-\i)A*[l/{\a^) + AW-\i)A*]-^y. (48) 

We demonstrate a procedure to compute C,{x) in (l48T l effi¬ 
ciently by using FFT. Let us rewrite (l48l l as 

C(x) = ^^-^(xjA*^ (49) 

where, y = [R + 1/(Acr^)]!:, 

R = AW'-i(x)A* 

Now partition x into G number of sub-vectors: 

- rzy»1 zir»T IT 

— ^2 d.Q\ 















where length of each Xi is Ni. Construct the matrix F by 
extracting first columns of the Fourier matrix F in (|4|i. We 
calculate ( |49] | by using the following steps. 

• At first we compute R = AW~^{x)A*. Since Wa[x) 
is a diagonal matrix, it follows using (fTTl) and (fOl) that 

G 

^ [0 diag(0Tcg)]FlF-i(ig)F*[0 diag(0Tcg)]T 


i=i 

= \Fw{xg)]k-t, (50) 

where w{xg) is a vector constructed from the diagonal 
components of W~^{xg). We compute Fw{xg) via FFT 
using N log 2 (A) floating point operations. Recall that the 
entries of the code matrix C are {+1,-1}. In addition 
0 is a row selector matrix. Hence, a multiplication 
by [0 diag( 0 '''cg)] does not require any floating point 
operation. In fact, constructing 

Rg := [0 diag(0Tc5)]FlF-i(5;3)F*[0 diag(0Tcg)]T 

needs to extract a small block of 'FW~^{xg)F*, and 
changing the signs of some entries. Finally, we compute 
R = which requires {G— l)(M + l)M/2 flops. 

Hence this step requires G{Alog 2 (A)+M(M+l)/ 2 } — 
M{M + l)/2 flops in total. 

• Calculate z by solving [R + I/(Acr^)] 2 : = y. By using 
Cholesky factorization, this M x M positive definite 
system of equations need 0{^M^) flops to compute z. 

m Compute A*z in parts, i.e. we partition A*z = 
[zi Z 2 ■ ■ ■ zg], and form 

Z 3 = F*[0 diag(0Tcg)]T2; 5=1,2,...G (51) 

by computing the IFFT of [0 diag( 0 '''Cg)]''' 2 :. Recall that 
forming [0 diag( 0 ‘'’Cg)]‘'’ 2 : does not require any addi¬ 
tional floating point operation. Hence this step requires 
0 (GAlog 2 (A)) flops. 

• Finally, as W~^{x) is diagonal, we need GNi multipli¬ 
cations to compute C(®)- 

Thus in total we need 0 ( 2 GAlog 2 (A) -f GM{M + l)/2 + 
GNi + 1/3M^(M — 1.5)) flops to compute C(i)- 


H. Thresholding the recovered signal from ISLO 

Let the final output obtained from ISLO is x. We denote 
V = X — X. The vector v can be viewed as the recovery error 
resulted due to noise. To perform a thresholding of x, we 
analyze the statistical property of v. Using (|4^ . we can write 

n ”1 


X = 


Wa{x) 

A(t 2 


A* A 


A*y. 


(52) 


Assuming ||u ||2 is small, the first order Taylor series expansion 
of Wc(x) around x is 

Wa{x) = Wa[x) - Wcr(a:)diag(-^)diag(u). (53) 


Now consider (l5^ . 

'W„{x) 


A*y = 


AtJ^ 


-f A*A 


Wa{x) 


Aa 2 


I-diag(^)-diag(|2)) +A*A 


j^W„ix)+A*A 
By ignoring second order terms in v we have 


( 54 ) 


W^ix) 

Aa 2 

= A*y- 


I-diag(—) -f A*A 


Aa 2 


Wa{x) + A*A 


X 


(55) 


For a small value of cr, we can neglect Wa{x)^, hence 


V = 

= De 




diag(^) -FA* A 


A*e 


(56) 


where we define ^ D = 

(l - diag(f^)) -F A*A A*. To compute D, we 
heed the value of x which is unknown in priori. Nevertheless, 
we can use an estimate of x to compute D. In this work, 
we use X as an estimate of x. Furthermore, computing D 
requires inverting a large size matrix. The computation task 
can be reduced significantly by applying matrix inversion 
lemma. Let us define P = 
verified that 


f I — diag(^) j. It can be 


D = P-^A* [I-F AP-^A*] \ (57) 

Let us partition the matrix D such that = 

[D(i) where each D« g Also 

partition v = ■ ■ ■ Vq]'^. Assume that e is complex 

Gaussian with zero mean and a covariance matrix Hence, 
the entries of D are independent of e. Then the variable 
11'*^*Hi = II II 2 has a generalized chi-square distribution 

of order M (assuming M < Ni) Il27ll . The procedure for 
computing the cumulative distribution function (CDF) of a 
variable having generalized chi-square distribution has been 
described in Il27l . Il28l . In this work, the CDF of ||wi ||2 for a 
threshold ti will be denoted by x(ri, where A^) is 

the vector containing the singular values of (Dd)[Dd)]*^. 

Partition x = [a:[ x^ ■ ■ ■ such that every Xi G C^L 
Define a set S' = {i : ||a;i ||2 ^ 0}. Then for a given i consider 
the two hypotheses: [TLq : i ^ S;7fi : i G S]. Note that under 
Ho, the distribution of ||Si||f is similar to the distribution of 
lirtilll. To perform hypothesis test on ||Si|| 2 , we need to select 
a threshold parameter r^. The procedure for selecting the value 
of Ti will be described next. The value of ||ii||| is checked 
against to take a decision between the two hypothesis: 



The threshold Ti is fixed to achieve a desired false alarm 
probability ip according to 


Ip = P{\\xi\\l > nlHo) 
= l-x(A,AW,cr^) 


(59) 
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There are total G number of sub-vectors i.e., The 

overall false alarm probability can be defined: 

P/a = 1 - (1 - iPf (60) 

To perform thresholding of x, we select a desired false 
alarm probability Pfa first. For the Pfa, we can calculate the 
threshold parameter for every ||aii||| by using ( |59] |. 

IV. Simulation Results 

A typical N = 1024 subcarrier OFDMA system, by 
following the WiMAX standards m, El, has been chosen for 
the simulation. In the system, the carrier frequency is 5.1 GHz, 
and the associated sampling interval is Tg = 89.28 ns. This 
corresponds to a subcarrier spacing of 10.94 kHz. Length of 
the cyclic prefix is 64 samples. Total M = 144 subcarriers are 
reserved for the initial ranging purpose, and the total number 
of available random codes in matrix C is 32, i.e. G = 32. 
The modulation pulse is a root-raised-cosine function with a 
roll-off 0.22 and duration lOTj. The channel impulse response 
has a maximum order P = 30, and the wireless cell radius 
is 2.5 km, hence D = 186. Similar to ||2l, 1^ . we assume 
that BS has an approximate knowledge about Pmax and we 
set Ni = Pmax + P in (flfil) . The RTs follows a mixed channel 
model specified by ITU IMT-2000 standards: Ped-A, Ped-B, 
and Veh-A. The RTs select the channel models with equal 
probability. The mobile speed varies in the interval [0,5] m/s 
for Ped-A, Ped-B channels, and [5,20] m/s for Veh-A. Since 
the ranging signal is used to measure the system performance, 
the signal to noise ratio is defined as SNR = lOlogj^g (^)’ 
where and are the variances of channel impulse response 
h and noise term e respectively. Four different algorithms are 
considered for performance comparison. The proposed algo¬ 
rithm will be called “Handover”. The other three algorithms 
are the SMUD ||2l, SRMD scheme discussed in a, and the 
MU-GLRT proposed in Il29ll . 

We start by finding a good choice of k to avoid unnecessary 
iterations in generating the rough estimate (to be used as the 
initial guess by ISLO) via the £i optimization. From a wide 
range of simulations with different number of IR users and 
SNR conditions we found that the performance of Handover 
remain almost same for k > 0.8. Hence we recommend setting 
K = 0.8, which is used in all the following cases. Figure 
[TJa) illustrates the code detection performance of the proposed 
algorithm for different values of false alarm probability Pfa 
(see dfiOl l). The performance is assessed in terms of success of 
code detection. Recall that the set of active IR code indices is 
£. Let £ be the set of code indices detected by an algorithm. 
The probability that £ = £, denoted by Pg, is used to quantify 
the merit of the algorithm We consider five different values 
of Pfa for code detection. As can be seen in Figure [11^), 
the Handover algorithm provides optimum performance for 
Pfa = le-4. Hence, we recommend setting Pfa = le-4. 

Figure [T]]b) shows the code detection performance by dif¬ 
ferent algorithms. The performance of MU-GLRT is worse 
for larger number of active ranging users compared to other 
three algorithms, whereas the Handover performs best. Note 
that at moderate SNR i.e., SNR= 10 dB, the Handover 


algorithm can recover 6 ranging users with high probability. 
The performance of SRMD is average compared to other 
algorithms. For instance, with 4 users and SNR= lOdB, 
the code detection probability of MU-GLRT, SRMD, SMUD 
and Handover are 0.36, 0.93, 0.91 and 0.98 respectively. The 
performance of MU-GLRT degrades rapidly with decreasing 
the SNR. Hence, we do not illustrate the result of MU-GLRT 
for lower SNR. With SNR= 3dB and 4 ranging users, the 
code detection probability of SRMD, SMUD and Handover 
are 0.88, 0.87 and 0.97 respectively. 

Figure |2] illustrates the accuracy of ranging parameter esti¬ 
mations by different algorithms. We do not compare the result 
with the MU-GLRT at low SNR as its performance is poor in 
the simulation environment. As can be seen in Figure |2l[ a), the 
MSB of power estimate increases with increasing the number 
of users. Note that with SNR= lOdB, the SMUD, MU-GLRT 
and Handover exhibit similar performance for small number 
of users (i.e, for total users 2 in Figure liJa)). However, their 
performance difference increases with increasing the number 
of users. The MSB of power estimate from SMUD for 2 
and 5 users (with SNR= 10 dB) are 0.0033 and 0.031 
respectively, whereas the MSB from Handover are 0.0034 
and 0.017 respectively. The MSB of the timing estimates for 
the considered ranging algorithms are shown in Bigure 12b). 
As can be seen, the Handover algorithm outperforms other 
algorithms with big margin. Bor example, with 4 ranging 
users and SNR= lOdB, the MSB for Handover is 4.9 which 
is 6.38 and 9.54 for SMUD and SRMD respectively. The 
MSB increases with increasing the number of users. Bor 
example, with 5 ranging users and SNR= 10 dB, the MSB 
of timing offset estimation by Handover, SMUD and SRMD 
are 6.24,9.132 and 12.27 respectively. We see that the SMUD 
performs better than SRMD at high SNR (10 dB), however 
SRMD outperforms SMUD at low SNR (i.e, 3 dB). We also 
compare computational complexity of the proposed algorithm 
with SMUD. The complexity of the SRMD algorithm has not 
been analysed in Q, hence we cannot incorporate the result 
in the figure. With SNR= 10 dB and total active users 6, 
the Handover and SMUD requires 2.58e7 and 1.64e7 flops 
in average respectively to resolve the IR request. 

V. Conclusion 

In this work, we explored a formulation of the OBDMA 
initial ranging parameter estimation problem in a sparse signal 
representation framework. We started with developing a math¬ 
ematical model that poses the ranging problem into a sparse 
signal recovery problem. An efficient procedure has been 
proposed that blends two different types of sparse recovery 
algorithms. The resulting algorithm exhibits efficient ranging 
parameter estimation performance. 
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